home *** CD-ROM | disk | FTP | other *** search
/ Sprite 1984 - 1993 / Sprite 1984 - 1993.iso / src / lib / m / RCS / support.c,v < prev   
Text File  |  1990-02-16  |  17KB  |  551 lines
  1. head     1.2;
  2. branch   ;
  3. access   ;
  4. symbols  ;
  5. locks    ; strict;
  6. comment  @ * @;
  7.  
  8.  
  9. 1.2
  10. date     90.02.16.12.52.34;  author douglis;  state Exp;
  11. branches ;
  12. next     1.1;
  13.  
  14. 1.1
  15. date     89.08.06.21.03.02;  author douglis;  state Exp;
  16. branches ;
  17. next     ;
  18.  
  19.  
  20. desc
  21. @original from BSD
  22. @
  23.  
  24.  
  25. 1.2
  26. log
  27. @fixes for ds3100
  28. @
  29. text
  30. @/*
  31.  * Copyright (c) 1985 Regents of the University of California.
  32.  * All rights reserved.
  33.  *
  34.  * Redistribution and use in source and binary forms are permitted
  35.  * provided that this notice is preserved and that due credit is given
  36.  * to the University of California at Berkeley. The name of the University
  37.  * may not be used to endorse or promote products derived from this
  38.  * software without specific prior written permission. This software
  39.  * is provided ``as is'' without express or implied warranty.
  40.  *
  41.  * All recipients should regard themselves as participants in an ongoing
  42.  * research project and hence should feel obligated to report their
  43.  * experiences (good or bad) with these elementary function codes, using
  44.  * the sendbug(8) program, to the authors.
  45.  */
  46.  
  47. #ifndef lint
  48. static char sccsid[] = "@@(#)support.c    5.2 (Berkeley) 4/29/88";
  49. #endif /* not lint */
  50.  
  51. /* 
  52.  * Some IEEE standard 754 recommended functions and remainder and sqrt for 
  53.  * supporting the C elementary functions.
  54.  ******************************************************************************
  55.  * WARNING:
  56.  *      These codes are developed (in double) to support the C elementary
  57.  * functions temporarily. They are not universal, and some of them are very
  58.  * slow (in particular, drem and sqrt is extremely inefficient). Each 
  59.  * computer system should have its implementation of these functions using 
  60.  * its own assembler.
  61.  ******************************************************************************
  62.  *
  63.  * IEEE 754 required operations:
  64.  *     drem(x,p) 
  65.  *              returns  x REM y  =  x - [x/y]*y , where [x/y] is the integer
  66.  *              nearest x/y; in half way case, choose the even one.
  67.  *     sqrt(x) 
  68.  *              returns the square root of x correctly rounded according to 
  69.  *        the rounding mod.
  70.  *
  71.  * IEEE 754 recommended functions:
  72.  * (a) copysign(x,y) 
  73.  *              returns x with the sign of y. 
  74.  * (b) scalb(x,N) 
  75.  *              returns  x * (2**N), for integer values N.
  76.  * (c) logb(x) 
  77.  *              returns the unbiased exponent of x, a signed integer in 
  78.  *              double precision, except that logb(0) is -INF, logb(INF) 
  79.  *              is +INF, and logb(NAN) is that NAN.
  80.  * (d) finite(x) 
  81.  *              returns the value TRUE if -INF < x < +INF and returns 
  82.  *              FALSE otherwise.
  83.  *
  84.  *
  85.  * CODED IN C BY K.C. NG, 11/25/84;
  86.  * REVISED BY K.C. NG on 1/22/85, 2/13/85, 3/24/85.
  87.  */
  88.  
  89.  
  90. #if defined(vax)||defined(tahoe)      /* VAX D format */
  91. #include <errno.h>
  92.     static unsigned short msign=0x7fff , mexp =0x7f80 ;
  93.     static short  prep1=57, gap=7, bias=129           ;   
  94.     static double novf=1.7E38, nunf=3.0E-39, zero=0.0 ;
  95. #else    /* defined(vax)||defined(tahoe) */
  96.     static unsigned short msign=0x7fff, mexp =0x7ff0  ;
  97.     static short prep1=54, gap=4, bias=1023           ;
  98.     static double novf=1.7E308, nunf=3.0E-308,zero=0.0;
  99. #endif    /* defined(vax)||defined(tahoe) */
  100.  
  101. double scalb(x,N)
  102. double x; int N;
  103. {
  104.         int k;
  105.         double scalb();
  106.  
  107. #if defined(national) || defined (ds3100)
  108.         unsigned short *px=(unsigned short *) &x + 3;
  109. #else    /* national  || ds3100 */
  110.         unsigned short *px=(unsigned short *) &x;
  111. #endif    /* national */
  112.  
  113.         if( x == zero )  return(x); 
  114.  
  115. #if defined(vax)||defined(tahoe)
  116.         if( (k= *px & mexp ) != ~msign ) {
  117.             if (N < -260)
  118.         return(nunf*nunf);
  119.         else if (N > 260) {
  120.         extern double infnan(),copysign();
  121.         return(copysign(infnan(ERANGE),x));
  122.         }
  123. #else    /* defined(vax)||defined(tahoe) */
  124.         if( (k= *px & mexp ) != mexp ) {
  125.             if( N<-2100) return(nunf*nunf); else if(N>2100) return(novf+novf);
  126.             if( k == 0 ) {
  127.                  x *= scalb(1.0,(int)prep1);  N -= prep1; return(scalb(x,N));}
  128. #endif    /* defined(vax)||defined(tahoe) */
  129.  
  130.             if((k = (k>>gap)+ N) > 0 )
  131.                 if( k < (mexp>>gap) ) *px = (*px&~mexp) | (k<<gap);
  132.                 else x=novf+novf;               /* overflow */
  133.             else
  134.                 if( k > -prep1 ) 
  135.                                         /* gradual underflow */
  136.                     {*px=(*px&~mexp)|(short)(1<<gap); x *= scalb(1.0,k-1);}
  137.                 else
  138.                 return(nunf*nunf);
  139.             }
  140.         return(x);
  141. }
  142.  
  143.  
  144. double copysign(x,y)
  145. double x,y;
  146. {
  147. #ifdef national
  148.         unsigned short  *px=(unsigned short *) &x+3,
  149.                         *py=(unsigned short *) &y+3;
  150. #else    /* national */
  151.         unsigned short  *px=(unsigned short *) &x,
  152.                         *py=(unsigned short *) &y;
  153. #endif    /* national */
  154.  
  155. #if defined(vax)||defined(tahoe)
  156.         if ( (*px & mexp) == 0 ) return(x);
  157. #endif    /* defined(vax)||defined(tahoe) */
  158.  
  159.         *px = ( *px & msign ) | ( *py & ~msign );
  160.         return(x);
  161. }
  162.  
  163. double logb(x)
  164. double x; 
  165. {
  166.  
  167. #ifdef national
  168.         short *px=(short *) &x+3, k;
  169. #else    /* national */
  170.         short *px=(short *) &x, k;
  171. #endif    /* national */
  172.  
  173. #if defined(vax)||defined(tahoe)
  174.         return (int)(((*px&mexp)>>gap)-bias);
  175. #else    /* defined(vax)||defined(tahoe) */
  176.         if( (k= *px & mexp ) != mexp )
  177.             if ( k != 0 )
  178.                 return ( (k>>gap) - bias );
  179.             else if( x != zero)
  180.                 return ( -1022.0 );
  181.             else        
  182.                 return(-(1.0/zero));    
  183.         else if(x != x)
  184.             return(x);
  185.         else
  186.             {*px &= msign; return(x);}
  187. #endif    /* defined(vax)||defined(tahoe) */
  188. }
  189.  
  190. finite(x)
  191. double x;    
  192. {
  193. #if defined(vax)||defined(tahoe)
  194.         return(1);
  195. #else    /* defined(vax)||defined(tahoe) */
  196. #ifdef national
  197.         return( (*((short *) &x+3 ) & mexp ) != mexp );
  198. #else    /* national */
  199.         return( (*((short *) &x ) & mexp ) != mexp );
  200. #endif    /* national */
  201. #endif    /* defined(vax)||defined(tahoe) */
  202. }
  203.  
  204. double drem(x,p)
  205. double x,p;
  206. {
  207.         short sign;
  208.         double hp,dp,tmp,drem(),scalb();
  209.         unsigned short  k; 
  210. #ifdef national
  211.         unsigned short
  212.               *px=(unsigned short *) &x  +3, 
  213.               *pp=(unsigned short *) &p  +3,
  214.               *pd=(unsigned short *) &dp +3,
  215.               *pt=(unsigned short *) &tmp+3;
  216. #else    /* national */
  217.         unsigned short
  218.               *px=(unsigned short *) &x  , 
  219.               *pp=(unsigned short *) &p  ,
  220.               *pd=(unsigned short *) &dp ,
  221.               *pt=(unsigned short *) &tmp;
  222. #endif    /* national */
  223.  
  224.         *pp &= msign ;
  225.  
  226. #if defined(vax)||defined(tahoe)
  227.         if( ( *px & mexp ) == ~msign )    /* is x a reserved operand? */
  228. #else    /* defined(vax)||defined(tahoe) */
  229.         if( ( *px & mexp ) == mexp )
  230. #endif    /* defined(vax)||defined(tahoe) */
  231.         return  (x-p)-(x-p);    /* create nan if x is inf */
  232.     if (p == zero) {
  233. #if defined(vax)||defined(tahoe)
  234.         extern double infnan();
  235.         return(infnan(EDOM));
  236. #else    /* defined(vax)||defined(tahoe) */
  237.         return zero/zero;
  238. #endif    /* defined(vax)||defined(tahoe) */
  239.     }
  240.  
  241. #if defined(vax)||defined(tahoe)
  242.         if( ( *pp & mexp ) == ~msign )    /* is p a reserved operand? */
  243. #else    /* defined(vax)||defined(tahoe) */
  244.         if( ( *pp & mexp ) == mexp )
  245. #endif    /* defined(vax)||defined(tahoe) */
  246.         { if (p != p) return p; else return x;}
  247.  
  248.         else  if ( ((*pp & mexp)>>gap) <= 1 ) 
  249.                 /* subnormal p, or almost subnormal p */
  250.             { double b; b=scalb(1.0,(int)prep1);
  251.               p *= b; x = drem(x,p); x *= b; return(drem(x,p)/b);}
  252.         else  if ( p >= novf/2)
  253.             { p /= 2 ; x /= 2; return(drem(x,p)*2);}
  254.         else 
  255.             {
  256.                 dp=p+p; hp=p/2;
  257.                 sign= *px & ~msign ;
  258.                 *px &= msign       ;
  259.                 while ( x > dp )
  260.                     {
  261.                         k=(*px & mexp) - (*pd & mexp) ;
  262.                         tmp = dp ;
  263.                         *pt += k ;
  264.  
  265. #if defined(vax)||defined(tahoe)
  266.                         if( x < tmp ) *pt -= 128 ;
  267. #else    /* defined(vax)||defined(tahoe) */
  268.                         if( x < tmp ) *pt -= 16 ;
  269. #endif    /* defined(vax)||defined(tahoe) */
  270.  
  271.                         x -= tmp ;
  272.                     }
  273.                 if ( x > hp )
  274.                     { x -= p ;  if ( x >= hp ) x -= p ; }
  275.  
  276. #if defined(vax)||defined(tahoe)
  277.         if (x)
  278. #endif    /* defined(vax)||defined(tahoe) */
  279.             *px ^= sign;
  280.                 return( x);
  281.  
  282.             }
  283. }
  284. double sqrt(x)
  285. double x;
  286. {
  287.         double q,s,b,r;
  288.         double logb(),scalb();
  289.         double t,zero=0.0;
  290.         int m,n,i,finite();
  291. #if defined(vax)||defined(tahoe)
  292.         int k=54;
  293. #else    /* defined(vax)||defined(tahoe) */
  294.         int k=51;
  295. #endif    /* defined(vax)||defined(tahoe) */
  296.  
  297.     /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
  298.         if(x!=x||x==zero) return(x);
  299.  
  300.     /* sqrt(negative) is invalid */
  301.         if(x<zero) {
  302. #if defined(vax)||defined(tahoe)
  303.         extern double infnan();
  304.         return (infnan(EDOM));    /* NaN */
  305. #else    /* defined(vax)||defined(tahoe) */
  306.         return(zero/zero);
  307. #endif    /* defined(vax)||defined(tahoe) */
  308.     }
  309.  
  310.     /* sqrt(INF) is INF */
  311.         if(!finite(x)) return(x);               
  312.  
  313.     /* scale x to [1,4) */
  314.         n=logb(x);
  315.         x=scalb(x,-n);
  316.         if((m=logb(x))!=0) x=scalb(x,-m);       /* subnormal number */
  317.         m += n; 
  318.         n = m/2;
  319.         if((n+n)!=m) {x *= 2; m -=1; n=m/2;}
  320.  
  321.     /* generate sqrt(x) bit by bit (accumulating in q) */
  322.             q=1.0; s=4.0; x -= 1.0; r=1;
  323.             for(i=1;i<=k;i++) {
  324.                 t=s+1; x *= 4; r /= 2;
  325.                 if(t<=x) {
  326.                     s=t+t+2, x -= t; q += r;}
  327.                 else
  328.                     s *= 2;
  329.                 }
  330.             
  331.     /* generate the last bit and determine the final rounding */
  332.             r/=2; x *= 4; 
  333.             if(x==zero) goto end; 100+r; /* trigger inexact flag */
  334.             if(s<x) {
  335.                 q+=r; x -=s; s += 2; s *= 2; x *= 4;
  336.                 t = (x-s)-5; 
  337.                 b=1.0+3*r/4; if(b==1.0) goto end; /* b==1 : Round-to-zero */
  338.                 b=1.0+r/4;   if(b>1.0) t=1;    /* b>1 : Round-to-(+INF) */
  339.                 if(t>=0) q+=r; }          /* else: Round-to-nearest */
  340.             else { 
  341.                 s *= 2; x *= 4; 
  342.                 t = (x-s)-1; 
  343.                 b=1.0+3*r/4; if(b==1.0) goto end;
  344.                 b=1.0+r/4;   if(b>1.0) t=1;
  345.                 if(t>=0) q+=r; }
  346.             
  347. end:        return(scalb(q,n));
  348. }
  349.  
  350. #if 0
  351. /* DREM(X,Y)
  352.  * RETURN X REM Y =X-N*Y, N=[X/Y] ROUNDED (ROUNDED TO EVEN IN THE HALF WAY CASE)
  353.  * DOUBLE PRECISION (VAX D format 56 bits, IEEE DOUBLE 53 BITS)
  354.  * INTENDED FOR ASSEMBLY LANGUAGE
  355.  * CODED IN C BY K.C. NG, 3/23/85, 4/8/85.
  356.  *
  357.  * Warning: this code should not get compiled in unless ALL of
  358.  * the following machine-dependent routines are supplied.
  359.  * 
  360.  * Required machine dependent functions (not on a VAX):
  361.  *     swapINX(i): save inexact flag and reset it to "i"
  362.  *     swapENI(e): save inexact enable and reset it to "e"
  363.  */
  364.  
  365. double drem(x,y)    
  366. double x,y;
  367. {
  368.  
  369. #ifdef national        /* order of words in floating point number */
  370.     static n0=3,n1=2,n2=1,n3=0;
  371. #else /* VAX, SUN, ZILOG, TAHOE */
  372.     static n0=0,n1=1,n2=2,n3=3;
  373. #endif
  374.  
  375.         static unsigned short mexp =0x7ff0, m25 =0x0190, m57 =0x0390;
  376.     static double zero=0.0;
  377.     double hy,y1,t,t1;
  378.     short k;
  379.     long n;
  380.     int i,e; 
  381.     unsigned short xexp,yexp, *px  =(unsigned short *) &x  , 
  382.                   nx,nf,      *py  =(unsigned short *) &y  ,
  383.                   sign,      *pt  =(unsigned short *) &t  ,
  384.                         *pt1 =(unsigned short *) &t1 ;
  385.  
  386.     xexp = px[n0] & mexp ;    /* exponent of x */
  387.     yexp = py[n0] & mexp ;    /* exponent of y */
  388.     sign = px[n0] &0x8000;    /* sign of x     */
  389.  
  390. /* return NaN if x is NaN, or y is NaN, or x is INF, or y is zero */
  391.     if(x!=x) return(x); if(y!=y) return(y);         /* x or y is NaN */
  392.     if( xexp == mexp )   return(zero/zero);      /* x is INF */
  393.     if(y==zero) return(y/y);
  394.  
  395. /* save the inexact flag and inexact enable in i and e respectively
  396.  * and reset them to zero
  397.  */
  398.     i=swapINX(0);    e=swapENI(0);    
  399.  
  400. /* subnormal number */
  401.     nx=0;
  402.     if(yexp==0) {t=1.0,pt[n0]+=m57; y*=t; nx=m57;}
  403.  
  404. /* if y is tiny (biased exponent <= 57), scale up y to y*2**57 */
  405.     if( yexp <= m57 ) {py[n0]+=m57; nx+=m57; yexp+=m57;}
  406.  
  407.     nf=nx;
  408.     py[n0] &= 0x7fff;    
  409.     px[n0] &= 0x7fff;
  410.  
  411. /* mask off the least significant 27 bits of y */
  412.     t=y; pt[n3]=0; pt[n2]&=0xf800; y1=t;
  413.  
  414. /* LOOP: argument reduction on x whenever x > y */
  415. loop:
  416.     while ( x > y )
  417.     {
  418.         t=y;
  419.         t1=y1;
  420.         xexp=px[n0]&mexp;      /* exponent of x */
  421.         k=xexp-yexp-m25;
  422.         if(k>0)     /* if x/y >= 2**26, scale up y so that x/y < 2**26 */
  423.         {pt[n0]+=k;pt1[n0]+=k;}
  424.         n=x/t; x=(x-n*t1)-n*(t-t1);
  425.     }    
  426.     /* end while (x > y) */
  427.  
  428.     if(nx!=0) {t=1.0; pt[n0]+=nx; x*=t; nx=0; goto loop;}
  429.  
  430. /* final adjustment */
  431.  
  432.     hy=y/2.0;
  433.     if(x>hy||((x==hy)&&n%2==1)) x-=y; 
  434.     px[n0] ^= sign;
  435.     if(nf!=0) { t=1.0; pt[n0]-=nf; x*=t;}
  436.  
  437. /* restore inexact flag and inexact enable */
  438.     swapINX(i); swapENI(e);    
  439.  
  440.     return(x);    
  441. }
  442. #endif
  443.  
  444. #if 0
  445. /* SQRT
  446.  * RETURN CORRECTLY ROUNDED (ACCORDING TO THE ROUNDING MODE) SQRT
  447.  * FOR IEEE DOUBLE PRECISION ONLY, INTENDED FOR ASSEMBLY LANGUAGE
  448.  * CODED IN C BY K.C. NG, 3/22/85.
  449.  *
  450.  * Warning: this code should not get compiled in unless ALL of
  451.  * the following machine-dependent routines are supplied.
  452.  * 
  453.  * Required machine dependent functions:
  454.  *     swapINX(i)  ...return the status of INEXACT flag and reset it to "i"
  455.  *     swapRM(r)   ...return the current Rounding Mode and reset it to "r"
  456.  *     swapENI(e)  ...return the status of inexact enable and reset it to "e"
  457.  *     addc(t)     ...perform t=t+1 regarding t as a 64 bit unsigned integer
  458.  *     subc(t)     ...perform t=t-1 regarding t as a 64 bit unsigned integer
  459.  */
  460.  
  461. static unsigned long table[] = {
  462. 0, 1204, 3062, 5746, 9193, 13348, 18162, 23592, 29598, 36145, 43202, 50740,
  463. 58733, 67158, 75992, 85215, 83599, 71378, 60428, 50647, 41945, 34246, 27478,
  464. 21581, 16499, 12183, 8588, 5674, 3403, 1742, 661, 130, };
  465.  
  466. double newsqrt(x)
  467. double x;
  468. {
  469.         double y,z,t,addc(),subc(),b54=134217728.*134217728.; /* b54=2**54 */
  470.         long mx,scalx,mexp=0x7ff00000;
  471.         int i,j,r,e,swapINX(),swapRM(),swapENI();       
  472.         unsigned long *py=(unsigned long *) &y   ,
  473.                       *pt=(unsigned long *) &t   ,
  474.                       *px=(unsigned long *) &x   ;
  475. #ifdef national         /* ordering of word in a floating point number */
  476.         int n0=1, n1=0; 
  477. #else
  478.         int n0=0, n1=1; 
  479. #endif
  480. /* Rounding Mode:  RN ...round-to-nearest 
  481.  *                 RZ ...round-towards 0
  482.  *                 RP ...round-towards +INF
  483.  *           RM ...round-towards -INF
  484.  */
  485.         int RN=0,RZ=1,RP=2,RM=3;/* machine dependent: work on a Zilog Z8070
  486.                                  * and a National 32081 & 16081
  487.                                  */
  488.  
  489. /* exceptions */
  490.     if(x!=x||x==0.0) return(x);  /* sqrt(NaN) is NaN, sqrt(+-0) = +-0 */
  491.     if(x<0) return((x-x)/(x-x)); /* sqrt(negative) is invalid */
  492.         if((mx=px[n0]&mexp)==mexp) return(x);  /* sqrt(+INF) is +INF */
  493.  
  494. /* save, reset, initialize */
  495.         e=swapENI(0);   /* ...save and reset the inexact enable */
  496.         i=swapINX(0);   /* ...save INEXACT flag */
  497.         r=swapRM(RN);   /* ...save and reset the Rounding Mode to RN */
  498.         scalx=0;
  499.  
  500. /* subnormal number, scale up x to x*2**54 */
  501.         if(mx==0) {x *= b54 ; scalx-=0x01b00000;}
  502.  
  503. /* scale x to avoid intermediate over/underflow:
  504.  * if (x > 2**512) x=x/2**512; if (x < 2**-512) x=x*2**512 */
  505.         if(mx>0x5ff00000) {px[n0] -= 0x20000000; scalx+= 0x10000000;}
  506.         if(mx<0x1ff00000) {px[n0] += 0x20000000; scalx-= 0x10000000;}
  507.  
  508. /* magic initial approximation to almost 8 sig. bits */
  509.         py[n0]=(px[n0]>>1)+0x1ff80000;
  510.         py[n0]=py[n0]-table[(py[n0]>>15)&31];
  511.  
  512. /* Heron's rule once with correction to improve y to almost 18 sig. bits */
  513.         t=x/y; y=y+t; py[n0]=py[n0]-0x00100006; py[n1]=0;
  514.  
  515. /* triple to almost 56 sig. bits; now y approx. sqrt(x) to within 1 ulp */
  516.         t=y*y; z=t;  pt[n0]+=0x00100000; t+=z; z=(x-z)*y; 
  517.         t=z/(t+x) ;  pt[n0]+=0x00100000; y+=t;
  518.  
  519. /* twiddle last bit to force y correctly rounded */ 
  520.         swapRM(RZ);     /* ...set Rounding Mode to round-toward-zero */
  521.         swapINX(0);     /* ...clear INEXACT flag */
  522.         swapENI(e);     /* ...restore inexact enable status */
  523.         t=x/y;          /* ...chopped quotient, possibly inexact */
  524.         j=swapINX(i);   /* ...read and restore inexact flag */
  525.         if(j==0) { if(t==y) goto end; else t=subc(t); }  /* ...t=t-ulp */
  526.         b54+0.1;        /* ..trigger inexact flag, sqrt(x) is inexact */
  527.         if(r==RN) t=addc(t);            /* ...t=t+ulp */
  528.         else if(r==RP) { t=addc(t);y=addc(y);}/* ...t=t+ulp;y=y+ulp; */
  529.         y=y+t;                          /* ...chopped sum */
  530.         py[n0]=py[n0]-0x00100000;       /* ...correctly rounded sqrt(x) */
  531. end:    py[n0]=py[n0]+scalx;            /* ...scale back y */
  532.         swapRM(r);                      /* ...restore Rounding Mode */
  533.         return(y);
  534. }
  535. #endif
  536. @
  537.  
  538.  
  539. 1.1
  540. log
  541. @Initial revision
  542. @
  543. text
  544. @d78 1
  545. a78 1
  546. #ifdef national
  547. d80 1
  548. a80 1
  549. #else    /* national */
  550. @
  551.